####################################
#Diff-in-diff with matching and perm test - function 
#Traffic and Political Attitudes
#Connor Jerzak & Brian Libgober
#see http://r.iq.harvard.edu/docs/zelig/3.4-8/Additional_Arguments_f3.html for matchit additional args 
####################################

dd_with_randomTest <- function(my_data, #data frame 
                       treatment_basis_name, #character string 
                       control_basis_name, #character vector 
                       outcome_time1_name, #character string 
                       outcome_time2_name, #character string 
                       matching_variables, #character vector
                       treatment_basis_treated_max, #numeric 
                       treatment_basis_control_min, #numeric 
                       control_basis_control_max, #numeric 
                       control_basis_treated_min, #numeric 
                       match_quietly=T, 
                       start_browser=F, 
                       max_perm_numb=10000000, 
                       my_method = "nearest", 
                       my_distance="mahalanobis", 
                       my_discard="none", 
                       m.order="largest", 
                       my_replace=T)
{#begin function brackets 
  require(MatchIt, quietly=T)
  require(ri, quietly=T)
  
  if(start_browser==T){browser()}
  
  #get treatment basis into memory 
  treatment_basis <- eval(parse(text=paste("my_data$", treatment_basis_name, "")))
  
  #get control basis into memory 
  control_basis <- eval(parse(text=paste("my_data$", control_basis_name, "")))
  
  #determine treated_indicator 
  my_data <- cbind(my_data, NA) 
  colnames(my_data)[ncol(my_data)] <- "treated_indicator"
  my_data$treated_indicator[treatment_basis <= treatment_basis_treated_max & 
                              control_basis >= control_basis_treated_min] <- 1
  my_data$treated_indicator[control_basis <= control_basis_control_max & 
                              treatment_basis >= treatment_basis_control_min] <- 0 
  summary(my_data$treated_indicator)
  
  #delete variable with is.na(treated_indicator)=T
  my_data <- my_data[is.na(my_data$treated_indicator)==F,]
  remove_columns <- colnames(my_data) %in% c(matching_variables, outcome_time1_name, outcome_time2_name, "treated_indicator")
  my_data <- my_data[,remove_columns]
  
  #create matching formula 
  matching_formula <- paste("treated_indicator ~", paste(matching_variables, collapse="+"))
  
  #create estimation_function 
  estimation_function <- function(match_input=my_data)
  { 
    #create matched dataset
    matched_results <- matchit(as.formula(matching_formula), 
                               data = match_input, 
                               replace=my_replace, 
                               m.order="smallest", 
                               method = my_method, 
                               distance=my_distance, 
                               discard= my_discard)
    
    #obtain matched dataset 
    matched_dataset <- match.data(matched_results)
    
    return(matched_dataset)
  }
  
  did_calc <- function(input_dataset, treated_index)
  { 
    mean_treat_time1 <- mean(input_dataset[treated_index==1,which(colnames(input_dataset)==outcome_time1_name)])
    mean_treat_time2 <- mean(input_dataset[treated_index==1,which(colnames(input_dataset)==outcome_time2_name)])
  
    mean_control_time1 <- mean(input_dataset[treated_index==0,which(colnames(input_dataset)==outcome_time1_name)])
    mean_control_time2 <- mean(input_dataset[treated_index==0,which(colnames(input_dataset)==outcome_time2_name)])
  
    diff_in_diff_estimate <- (mean_treat_time2 - mean_treat_time1) - (mean_control_time2 - mean_control_time1)
    return(diff_in_diff_estimate)
   }
  
  ####################
  #calculate base result 
  ####################
  base_matched_dataset <- estimation_function(match_input=my_data)
  base_result <- did_calc(input_dataset=base_matched_dataset, 
                                          treated_index=base_matched_dataset$treated_indicator)
  
  ############
  #permutation test 
  ############
  permutation_matrix <- genperms(base_matched_dataset$treated_indicator, blockvar = NULL, clustvar = NULL, maxiter = max_perm_numb)

  perm_results <- apply(permutation_matrix, 2, function(x) did_calc(input_dataset=base_matched_dataset, treated_index=x))

  #find the p-value 
  p_value <- mean(abs(perm_results) >= base_result)
    
  return_list <- list(base_result=base_result,
                      p_value=p_value)

  return(return_list)
}#end function barkets 
